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Abstract 

We model the dynamics of the leaky integrate-fire neuron under periodic 
stimulation as a Markov process with respect to the stimulus phase. This 
avoids the unrealistic assumption of a stimulus reset after each spike made 
in earlier work and thus solves the long-standing reset problem. The neuron 
exhibits stochastic resonance, both with respect to input noise intensity and 
stimulus frequency. The latter resonance arises by matching the stimulus 
frequency to the refractory time of the neuron. The Markov approach can 
be generalized to other periodically driven stochastic processes containing a 
reset mechanism. 
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I. INTRODUCTION 



Periodically modulated stochastic processes have been studied intensely over the last two 
decades under the paradigm of stochastic resonance: the transduction of signals is optimal 
in the presence of a particular amount of noise. First suggested to explain the periodicity of 
ice-ages |3J, stochastic resonance has since been demonstrated in a wide range of experiments 
and the underlying mechanisms are well understood. A recent review of the field is given 
in@. 

The concept of stochastic resonance has met with particular attention in the neuro- 
sciences • The brain achieves an enormous signal processing performance in the presence 
of noise from a wide range of sources, ranging from stochastic membrane channel openings 
on a molecular level, via highly irregular firing patterns of individual neurons to distracting 
stimuli in perception. The improvement of signal transduction on all of these levels has now 
been demonstrated experimentally P~|TH- Recently, the first direct evidence for the behav- 
ioral relevance of stochastic resonance has been reported underlining the importance 
of stochastic resonance in neurobiology. 

In short, neurons are threshold devices that receive an input I(t) which charges the 
membrane of the neuron like a leaky capacitor. When the potential v(t) across the membrane 
reaches a threshold B, a spike is fired: the membrane potential makes a brief but strong 
excursion (duration rs 2ms, amplitude ss 100m V). This spike is transmitted as output to 
other neurons. After the spike, the membrane potential is reset to a resting value Vq, some 
30mV below the threshold fll2"f . As the shape of the spikes is stereotypical, information is 
only conveyed by the spike times. 



This has led to the leaky integrate-fire model of neuronal dynamics . In between two 
spikes, the membrane potential is governed by 



T, 



l v(t) = -v(t) + I(t) + ((t). (1) 



Here, r m is the time-constant of the membrane, which represents the internal time-scale of the 
neuron and ((t) is an as yet undefined noise process, comprising, e.g., stochastic membrane 
potential fluctuations and irregular input to the neuron from sources uncorrelated to I(t). 
As the potential reaches the threshold, a spike is recorded and the potential is reset to 
v(t) = v o instantaneously. 

For Gaussian white noise ((t) the evolution of the membrane potential v (t) from reset 
potential to threshold is equivalent to an Ornstein-Uhlenbeck process with drift I(t) and an 
absorbing boundary at v — ©. The output of the neuron is modeled as a sequence of delta 
pulses fit) = J2kd(t ~ tk) at the times of threshold crossings {tk} = {t\v(t) = 0} (spike 
train). The spike train is a stochastic point process, specified entirely by the spike times 
{t k }. 

This biologically most interesting stochastic process has so far escaped a rigorous anal- 
ysis, in spite of several partially successful attempts [I^,I5|. For a list of open issues see 



Sec. V.C.4 of the review by Gammaitoni et al. 0. This is in marked contrast to the treatment 
of mathematically more accessible, but biologically less plausible models, such as bistable 
dynamic systems [|16HT9[1 and threshold devices without reset p0|-f22]|, in which stochastic 



resonance has been well established. 
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The essential difficulty arises from the reset after each spike: there is no well-defined 
membrane potential distribution for asymptotic times, as used in the case of reset-free 
threshold detectors. Instead we have to analyze each inter-spike-interval separately and then 
put these pieces together to obtain the spike train as a whole. To facilitate this, past work has 
assumed that the durations of all inter-spike-intervals (rfe = tu — tk-i) were identically and 
independently distributed (i.i.d.), i.e. that the spike train is a stationary renewal process |23 



But in the presence of time-dependent input I(t), this would require identical input within 
each inter-spike-interval (ISI). This is the much criticized reset assumption: If the model 
neuron were to describe a neuron in the auditory nerve while you are listening to a music 
tape, the reset assumption requires that upon the firing of each spike the tape should be 
rewound to exactly the position it had at the time of the last spike! 

In this work we show how to analyze the response of the leaky integrate-fire neuron to 
periodic stimuli without undue assumptions. The distribution of the length of individual 
inter-spike intervals is computed numerically ||15|| , and spike trains are then assembled as 
Markov chains from these intervals. We obtain probability distributions for the length 
of inter-spike intervals and the stimulus phases at which spikes occur. These distributions 
should be directly comparable to experiments employing sustained stimulation with periodic 
signals. The signal processing performance of the neuron is judged by the signal-to-noise 
ratio (SNR) of the output spike train. The SNR is maximal at an optimal noise amplitude 
for fixed stimulus frequency and at a resonance frequency for fixed noise amplitude. The 
latter resonance is a consequence of a time-scale matching between stimulus and membrane 
time-constant. All computations are verified by simulations. 

In Sec. [TI|, we show how to exploit the Markov property of the integrate-fire neuron to 
determine its response to sinusoidal input I(t). The performance of the model neuron as a 
signal processing device is investigated in Sec. |TTT|. The results are discussed in Sec. [TV| 



II. MARKOV ANALYSIS 

For an input current consisting of a constant offset and a sinusoidal component, and 
Gaussian white noise the Langevin equation (|I|) reads 

v(t) = -v(t) +/J + q cos(fit + O ) + VD£{t) , (2) 

where time and potential have been scaled to their respective natural units r m and O; the 
reset potential is set to v = 0. The input is characterized by the DC offset //, stimulus 
amplitude q, frequency f2 and initial phase <p . The noise has amplitude \f~D and auto- 
correlation (£(£)£(£')) = — ^')- m the remainder of this article, we will investigate this 
model. For a derivation of the type of input current used here from more elementary models, 
sec 



23| 



In the absence of noise (D — 0), spikes will only be generated if 

= lim v(t) = fi+ q > 1 . 

Therefore, we classify stimuli as sub-threshold if < 1 and as supra-threshold otherwise. 



In this work, we will focus on the biologically more interesting sub-threshold regime [25 
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The methods presented here are applicable independent of the choice of stimulus parameters. 
We only require the presence of noise, i.e. D > 0. 

Suppose that an initial spike has occured at time t = 0, corresponding to stimulus 
phase 0o- The next spike follows at time t\ = inf{t > to\v(t) > 1} and stimulus phase 
0i = (flti + 0o ) m °d 27r, whence cos(f2(t — ti) + 0i) = cos(f2t + O ) for t > t\. This suggests 
to re-write Eq. @ in terms of the time t' that has passed since the most recent spike at 
phase 0. Thus, given this phase, the potential evolves from v(t' = O|0) = vq = until the 
next spike according to 

v(t'\<p) = -v(t'\<p) +n + q cos{Qt' + 0) + y/D£(f) . (3) 

The next spike is fired after an interval r, as soon as the threshold condition is met 

r = inf{t' > 0\v(t'\(j)) > 1} . (4) 

The inter-spike intervals are connected by the iteration equations 

fc = (Qr k + fc _i) mod 27r , t k = t k -i + r k , (5) 

leading to the output spike train 

CO oo j 

/w=E^-*i) = E*(*-E^)- (6) 

j=0 j=0 k=l 

The reset of the membrane potential to vq = after each spike completely erases the 
memory of the neuron. The subsequent behavior of the neuron therefore depends on its past 
only through the absolute time of the spike t k , i.e. the spike train is a Markov process. 

We have thus split the task of solving the dynamics of the integrate-fire neuron into 
two parts. We will first solve the first-passage-time problem posed by Eqs. (Q) and (P for 
a given phase of the last spike, before assembling the spike train from the inter-spike 
intervals according to Eqs. @ and @. 



A. Conditional ISI distribution 

The first-passage-time problem for the membrane potential posed by Eqs. (|3], |j) yields 
the distribution p(r|0) of the inter-spike-interval lengths r for a given stimulus phase at 
the beginning of the interval (conditional ISI distribution). To the best of our knowledge, 
no analytic solution is known for this seemingly simple first-passage time problem of the 
Ornstein-Uhlenbeck process. The approximations suggested in |T4j are valid in a restricted 



parameter range only — low stimulus frequencies in particular — and appear to yield qualita- 
tive rather than quantitative agreement with simulations. 

We employ here a numerical method to compute the inter-spike-interval distributions. 
The method is discussed in detail in ||15|1 , and we only sketch it here. In the absence of an 
absorbing threshold the probability V(w, t'\u, s'; 0) that the membrane potential is w at time 
if if it was u at time s' < t' is a Gaussian distribution. The mean is given by the solution 
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at time t' of Eq. (|3]) for the noise- free case (D = 0) with initial condition v(s'\(p) = u, while 
the variance is cr 2 (t') = y(l — e~ 2<yt ~ s Then, the inter-spike-interval distribution is given 
by the integral equation []56| 

P(l,t'|0,0)= f* V(l,f\l,T) P (T\(i>)dT . (7) 

Jo 

This equation is solved for p using standard techniques J^fl. Source code is available on 
request. 

As shown in Fig. [TJ, the conditional inter-spike-interval distributions p(r\(p) may depend 
strongly on <p. First, they contain a series of exponentially decaying peaks that are separated 
by the stimulus period T = 2n/Q. These peaks represent spikes that are well phase- locked to 
the stimulus and we will refer to them as periodic peaks. An additional peak appears at short 
intervals r for certain phases (p. This peak reflects the rise time of the membrane potential 
towards threshold. Its location is not related to the stimulus period T, but reflects the 
intrinsic time-scale of the neuron and defines its refractory time, i.e. the minimum interval 
between two spikes. Thus, we will refer to this peak as the refractory peak. It corresponds 
to two or more spikes fired in rapid succession within a single stimulus period (a burst). 
There is thus a qualitative dependence of the distributions p on the phase that can lead 
to interesting consequences for the firing behavior of the neuron. 

The relation between periodic and refractory peaks depends on the stimulus parameters, 
particularly on the frequency and the noise amplitude. We will discuss this relationship in 
Sec. |TC|. ' 

B. Markov process in phase 

Let us now turn to the problem of assembling spike trains from inter-spike-intervals 
according to Eqs. ([|, ^|). The length of an interval following a spike at time t and stimulus 
phase (ft = [fit + O ] mod 27r is distributed according to p(r\(p). Therefore, the probability 
that the next spike will occur at phase ip is given by 

f°° dr 
T(ip\(f>)= / p{r\(j>)8{^- + mod27r) — . (8) 

Jo ' ' 

We will call T (ip | <fi) the transition probability of the spike phase. We will now consider the 
Markov process of the spike phases <ftk instead of the Markov process made up of the spike 
times tfc. 

If we define the spike phase distribution x (0) as the probability (across an ensemble 
of neurons or repetitions of an experiment) that the k th spike in a train will be fired at 
stimulus phase 4>, then this probability will evolve according to 

X (fc+1) W= / 2 V(^|0) X ( fc )(0)#. (9) 
Jo 

As the neuron fires repetitively while driven by a stationary periodic stimulus, the spike train 
emitted by the neuron will approach a stationary Markov process with phase distribution 
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k- 



2tt 



X (s) (V) = lim X (k) W= T(ip\ <P) X (s) W dcj> . (10) 



The stationary phase distribution x C0) * s the eigenfunction to eigenvalue 1 of the kernel 
T {ip | 0) , and is guaranteed to exist because this kernel is a conditional probability distri- 
bution p8[ . Any initial phase distribution will converge to the unique stationary solution 



provided that T (-0 | </>) > everywhere . That the latter condition holds in the presence 
of noise can be seen as follows. For sub-threshold stimuli, noise may drive the potential 
across the firing threshold at any time r > in principle, yielding a possibly tiny, but non- 
zero probability of spikes at any phase. The same argument holds true for supra-threshold 
stimuli, where noise may keep potential below threshold up to any time. In the absence of 
noise, neither convergence nor uniqueness are assured. 

To facilitate numerical treatment, we discretize the phase. Since the conditional inter- 
spike- interval distributions p(t|0) are smooth in both time and phase due to the presence of 
noise in the input, this discretization will introduce only minor numerical errors. It is largely 



equivalent to applying numerical methods to solve the kernel eigenvalue problem ||28|| . Using 
L bins of width Aip (Aip = 2tt/L) we obtain the spike phase distribution vector 

X= (Xo,Xi,- ■ ■ ,XL-i) tV , Xj = x(V0# , (ii) 

and the phase transition matrix T with elements 

Hj+i)AV> 

T lk = / T(^|A;AV0#, j,k = 0,...,L-l. (12) 

The evolution equation @ simplifies from convolution to matrix-vector multiplication 

x (k+i) = T . x (k) > (13) 

and the stationary distribution x 1S the eigenvector to eigenvalue 1 of the matrix T. We 
have thus reduced the Markov process to a Markov chain. 

In practice, we obtain the transition matrix T by numerically evaluating equations (§) 
and (0), with p(r\<p) from Eq. ([?]). The stationary distribution is then found using standard 
eigenvector routines. For all data shown here, we used the discretization L = 72, Atp = 
7r/36 = 5°. In figures of transition matrices and phase distributions the axis will run from 
—ii to 7i as this renders structures more clearly. 

An example for the phase evolution of an initially uniform distribution towards the 
stationary state under the influence of a transition matrix T is given in Fig. |. To "read" 
the transition matrix, note that the matrix columns correspond to the phase (f>k of the spike 
preceding the interval, the rows to the phase (fik+i of the spike terminating it. The phase 
axes run from — ti to 7r from bottom to top in phase distribution vectors x an d the rows of 
the transition matrix T, and from right to left across the columns of T. Thus, the horizontal 
bar in the transition matrix shown in Fig. ^| indicates that for most values of <pk the next 
spike will occur around (fik+i ~ — vr/6. This bar corresponds to the periodic peaks of the ISI 
distributions. For — 7r/4 < (f>k ^ tt/6, the matrix is dominated by a "finger", running parallel 
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to the matrix diagonal. Within this range of phases, a spike will be followed by another 
spike at a slightly later phase, as shown in Fig. ^b. Figuratively speaking, the neuron fires 
a burst of spikes, but there is always a chance that two subsequent spikes will be one or 
more stimulus periods apart, even though they are close in phase: in the Markov chain 
description, all information about actual interval lengths is lost. The finger results from the 
refractory peak of the ISI distributions. 

Figure |3| shows the dependence of transition matrix and stationary phase distribution 
on the noise amplitude for slow stimuli (T > 10). For low noise, the transition matrix 
is dominated by the horizontal bar, which intersects with the matrix diagonal, indicating 
a stochastic fixed point. This results in a sharply peaked spike phase distribution. At 
intermediate noise, the finger is more pronounced, while the bar barely touches the matrix 
diagonal, leading to a stochastic limit cycle with two preferred phases: the neuron often 
fires bursts of two successive spikes. At high noise, the finger stretches all along the matrix 
diagonal, while the horizontal bar has disappeared altogether. The neuron fires rapidly, but 
largely uncorrelated with the stimulus and the phase distribution is virtually flat. 

This means that for very low noise the spike train of the neuron is nearly a stationary re- 
newal process with inter-spike-intervals i.i.d. according to p{r\ip*). Here tp* is the location of 
the maximum of the stationary phase distribution, which depends not only on the stimulus 
parameters, but also on the noise amplitude. For high noise amplitudes, the response of the 
neuron is largely independent of the stimulus, and may thus be described by a stationary 
renewal process as well — the ISIs reduce to the refractory peak. But at intermediate noise 
levels — i.e. those essential to the observation of stochastic resonance — the stationary phase 
distribution may be multimodal. Thus the correlations between the phases of subsequent 
spikes have to be taken into account using the Markov ansatz. Multimodal phase distribu- 
tions as discussed here are not just hypothetical: they have been observed in sensory neurons 
of goldfish upon stimulation with sinusoidal water waves . 

For fast stimuli (T < 10), the stationary phase distribution smears out much more along 
the phase axis, and does not show multimodality, because the refractory time of the neuron 
becomes comparable to the stimulus period and bursting is no longer possible, see Fig. |j. At 
low to intermediate noise, the distribution is too wide to be replaced by its mode as in the 
renewal ansatz, but still sufficiently narrow to provide for a response that is well phase-locked 
to the stimulus. Therefore, the Markov approach is essential for high frequency stimuli as 
well. 



C. Stationary ISI distribution 

Once the stationary phase distribution is known, the inter-spike-interval distribution of 
the stationary firing process is obtained by averaging the conditional ISI distributions over 
phase 

r-2-K 

p(r)= / p(r|V)x (s) W#- (14) 
Jo 

The average interval length thus is 
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(r) = / rp{r)dr . (15) 



o 

p(r) is the inter-spike-interval distribution that we expect to find in experiments with tonic 
stimulation. In contrast to a stationary renewal process, this averaged ISI distribution does 
not contain a full description of the spike train. 

Typical ISI distributions p(r) are given in Figs. [5] and [5] for the same parameters as used 
in Figs. 0, H, respectively. For low noise, they contain only periodic peaks, located precisely 
at integer multiples of the stimulus period T: the neuron can only fire in a small time window 
within each period, and several periods may be skipped in between spikes. This indicates a 
firing pattern that is well phase-locked to the stimulus. ISI distributions with comparable 
structure have been found in neurons of the auditory system in different species pT| . |32f . For 
high noise, the ISI distributions reduce to the refractory peak, i.e. a largely random firing 
pattern. 

For intermediate noise, the ISI distributions depend strongly on the stimulus frequency. 
For high frequency (Fig. |B|), we find merely a superposition of periodic and refractory peaks: 
spikes preferentially occur at intervals that are multiples of the stimulus period, but this 
phase-locking is weak. This is very different for slow stimuli (Fig. ||), where the refractory 
peak is clearly separated from a wide peak at r = T = 40, the latter exposing some sub- 
structure. This can be understood as follows. The maximum of p(r) at r = T corresponds 
to two spikes fired each at the optimal phase in two subsequent periods. In contrast, if a 
period that contained a burst of two spikes is followed by another period containing a burst, 
then typically the first spike will be slightly earlier than the optimal phase, the second one 
a bit later. Thus, the interval between the second spike of the first burst and the first spike 
of the second burst is shorter than the stimulus period, leading to the side-peak at r ~ 35. 
The bursts themselves give rise to the refractory peak. This again indicates that the spike 
train is not a stationary renewal process. 

Along with results obtained using the Markov chain approach, Figs. |^-^ display phase and 
ISI distributions obtained from simulated trains of 20,000 spikes. The agreement between 
Markov model and simulation is excellent. Source code for the simulation based on [ 33 1 is 
available on request. 



III. STOCHASTIC RESONANCE 

To assess the performance of the integrate-fire neuron as a signal processing device, we 
evaluate the signal-to-noise ratio (SNR) of the spike train generated in response to periodic 
input. In doing so, one should keep in mind the purpose of the output spike train. It has to 
convey information to other neurons in the brain within a certain time window, as the brain 
has to respond quickly to stimuli. Therefore, the relevant quantity is the signal-to-noise ratio 



that can be achieved by measuring the spike train over a finite observation time T Q ||34j| . 
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A. Signal-to-noise ratio 



The one-sided power spectral density of a stationary spike train f(t) [as defined in Eq. (^j)] 
over a time interval T Q is |35| 

° f(t)e iu3t dt 



(16) 



S' To {uj)-- 




The average is to be taken over the ensemble of all spike trains, that is, over the set of 
all conditional ISI distributions and their (j — fc)-fold convolutions. This problem appears 
intractable. 

The situation is greatly simplified if u is the stimulus frequency Q or one of its harmonics. 
Expressing the spike times as tj = (rrij + p-)T, Eq. (|16"D for uo = nVt simplifies to 



S' To (nVL) 




(17) 



where n, rrij are integers, ipj G [0, 27r), and T = 2it/Q is the stimulus period. In the obser- 
vation period T a , on average T j (r) spikes will occur, regardless of the detailed structure of 
the spike train. We therefore fix the upper limit of the summation at N a = \T Q / (r)J , where 
(r) is the average interval length from Eq. flT5p and [x\ is the largest integer not exceeding 
x. This yields as an approximation 



S' To (ntt) w S To {ntt) 



7lN (t) 




>wj-</>t) 



(18) 



The task of computing an expectation with respect to all possible spike trains is now 
reduced to that of averaging over all possible sequences of spike phases. Their distribution 
and correlations are completely characterized by the transition matrix T, permitting for 
evaluation of Eq. ([18]) in closed form. The actual calculation is straightforward albeit lengthy 
algebra and is provided in the appendix. The final result may be written as 



s To (nn) 



tt(t) 



1 + A(n,N ) + {N -l)B(n) 



(19) 



where the functions A(n,N ) and B(n) are given in the appendix. Note that A(n, N a ) is 
bounded as iV — > oo. For a Poissonian spike train, both A and B are identically zero, 
yielding a white power spectrum p| . 

At first, it might seem surprising that the spectrum contains a term, (N — l)B(n), 
that scales linearly with the number of spikes in the train. This is a consequence of the 
periodic component of the spike train introduced by the driving stimulus, leading to a mixed 



spectrum consisting of a continuous background and a discrete spectrum of harmonics [35 
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For infinite observation time, i.e. N a — > oo, this gives rise to the terms ~ 8{u — nQ) in the 
power spectrum. 

A typical power spectrum is shown in Fig. |7|, indicating close agreement of Eq. ([19]) 
with results obtained by numerical Fourier transformation of simulated spike trains. The 
approximation made in fixing the summation limit in Eq. (|T8| ) is therefore well justified. 
The dip in the noise background of the spectrum at low frequencies is a consequence of the 
refractory period of the neuron, while the weak hump at u ~ 1 indicates the presence of 
bursts ||36|| . Spectra consisting only of this background have been found in neurons of higher 
cortical areas of monkeys in the absence of periodic input 

Since the power spectral density can only be evaluated in closed form at multiples of 
the stimulus frequency, we approximate the noise background as Poissonian white noise 
Sp = (ti (t)) _1 of a spike train of equal intensity p4"| . The signal-to-noise ratio obtainable 
from the spike train within the observation time T a is therefore given by 



SNR 



To 



Sp 



1+4(1, 



(r)J 



+ 



Ik 
(r) 



DB(1) 



(20) 



The signal-to-noise ratio for three different stimulus frequencies is shown in Fig. ^ vs. the 
noise amplitude, again in excellent agreement with simulation results. Stochastic resonance 
(SR) is clearly present at all frequencies, as the SNR attains its maximum for an intermediate 
noise level. The striking new feature is that the overall maximum in the SNR is reached 
at an intermediate frequency Q r ~ 7r/3, which we thus call the resonance frequency. The 
same qualitative dependence of the SNR on noise amplitude and stimulus frequencies is 
observed over a wide range of stimulus parameters, including weakly supra-threshold cases 
(0.4 < /I < 1, 0.4 < q/{\ - fj) < 1.2; data not shown). 

Note that the stochastic resonance reported in an earlier paper |15] is an artifact of the 
renewal ansatz employed in that work. There, the stimulus phase is reset to an arbitrarily 
chosen value 0o after each spike, and the signal-to-noise ratio is computed for an infinite 
observation time. The SNR is maximized for that noise level at which the periodic peaks 
of the ISI distribution p(r |0o) are centered about the multiples of the stimulus period T. 
But if, for low noise, one uses for each noise level D a different <po(D), namely the mode 
of the stationary phase distribution as discussed in Sec. 11 B , the periodic peaks are at 
multiples of T for all noise intensities, whence the SNR does not drop off for D —>■ and 
no resonance occurs (data not shown). This observation underlines the importance of the 
Markov approach. 



B. Time-scale matching 

In contrast to stochastic resonance in dynamical systems, SR with respect to the noise 
amplitude is not induced by the matching of time-scales in threshold systems, but results 
from stochastic linearization of the response function of the neuron In contrast, 

the additional resonance along the frequency axis arises in the integrate-fire neuron as a 
consequence of matching the stimulus period to the intrinsic time scale of the neuron in 
an appropriate manner. This is demonstrated in Fig. |9|. For a stimulus at the resonance 
frequency Q r , the peak at r = T in the stationary ISI distribution can "grow" in place as 
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noise is increased, without being disturbed by the refractory peak. Indeed, the latter arises 
at the location of the first periodic peak and shifts away from r = T only for very large noise. 
In this way, the firing rate of the neuron can be increased without loosing the phase-locking 
to the stimulus. Compare this to the cases of lower (Fig. ||) and higher (Fig. ^) frequencies: 
in both cases, high firing rates can only be achieved by raising the noise amplitude to a 
point where the refractory peak has either replaced (Q < Q r ) or smeared out (Q > fl r ) the 
periodic peaks, resulting in a firing pattern poorly phase-locked to the stimulus. 

This competition of precision and intensity is demonstrated by a phenomenological ansatz 
for the SNR. A measure of phase-locking between stimulus and response is the vector 
strength C s = (e 1 ^) , where ip are the spike phases p?[ . C s — 1 indicates perfect and 
C s = no locking. If the neuron attempts to measure the degree of phase-locking from a 
train of N = T Q j (r) spikes, the quality of measurement will be ~ y/N. Thus, we expect 
that the signal-to-noise ratio will roughly given by 



SNR phen « C S VN = cj-^. (21) 

Figure [TOJ demonstrates that this simple model describes the behavior of the SNR well. In 
particular, the two-fold stochastic resonance is reproduced. 

In short, to elicit a strong output signal from the model neuron, a sufficient input noise 
level is required. But this comes at a cost, as the quality of the output, i.e. the precision of 
the phase locking, deteriorates as noise is added. The maximum SNR represents the optimal 
compromise between signal strength and quality. 



IV. DISCUSSION 

In this paper, we have shown that the periodically driven integrate-fire neuron can be 
analyzed in the framework of a Markov process. This avoids the unrealistic assumption of a 
stimulus reset after each spike, the most serious shortcoming of previous work ||l4| , |l5| , and 



this answers question (1) raised by Gammaitoni et al. in Sec. V.4.C of their review 0. Their 
second questions concerns the fact that the neural membrane is a rectifier: even a strong 
negative input current will not lower the membrane potential more than a few millivolts 
below the reset potential v Q . This would indeed be a problem if the DC offset /i of the 
input were much smaller than the amplitude q of the AC stimulus. Preliminary evidence 
suggests that the best fit of inter-spike-interval distributions generated by the model with 
experimental data from the cat's auditory system is obtained for sub-threshold stimuli 



with /i > g. In this regime, the membrane potential is quickly raised to vq + fi and then 
oscillates around this level, unaffected by rectification. Finally, Gammaitoni and co-authors 



question the validity of the approximations used to compute the ISI distributions in 
This matter is avoided here by numerically computing these distributions. A study of the 
validity of approximate closed-form ISI distributions will be given elsewhere [50]. 



The Markov formalism presented in this paper is applicable to any periodically driven 
stochastic process with a reset. The only required ingredients are the conditional first- 
passage-time distributions p(r\(f>) and the iteration equations @. The generalization to 
more complex stimuli, e.g. including amplitude modulation, is straightforward. 
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With the Markov machinery at hand, we have demonstrated that the signal-to-noise ratio 
of the output of the neuron is maximized at an optimal noise amplitude for fixed frequency 
and at a resonance frequency for fixed noise intensity. Stochastic resonance with respect to 
the stimulus frequency, termed bona fide stochastic resonance, has been described in bistable 
systems before [[Il],[|IJ. Therefore, our findings for a non-dynamical threshold neuron extend 
the universality of stochastic resonance to the case of bona fide SR. Recent criticism |^3| of 
the original definition of bona fide SR, based on residence time distributions, does not apply 
to our study. 

Neurons in the auditory system can phase lock to acoustic stimuli with high acuity and 



utilize this for the precise localization of sound sources [[44] . Our results show that strong 
signals that are well phase locked to a stimulus may be achieved in spite of the noise ubiq- 
uitous in the neural system. Stochastic resonance might therefore be one of the underlying 
mechanisms of stereo hearing. First qualitative comparisons indicate good agreement be- 
tween response properties of the integrate-fire neuron and of auditory neurons. An intriguing 
question in this respect is the relevance of the bona fide SR to the neural system. It may 
serve to tune neurons as bandpass filters of a special kind: only stimuli in a certain frequency 
window will be transmitted with high intensity and precise phase locking. A detailed study 
will be the topic of a future publication. 
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APPENDIX A: COMPUTING THE POWER SPECTRAL DENSITY 

To prove Eq. (|I9|) , i.e. 



X 1 \j,k=l 



M 

in(ipj-il> k ) 



1 + A(n, M) + (M-l)S(n) 



7r(T) 

we split the double sum into the diagonal and off-diagonal terms 

S To (nQ) = [1 + h M (nQ) + h* M (nn)\ , (Al) 

7T(t) 

1 M M-k 

h M (nQ) = — E <e m{ ^-^) , (A2) 

k=l 3=1 
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the asterisk denoting complex conjugation and M = \T d /t\. 

Since we are considering a stationary Markov process, all ipk are identically distributed 
according to x > while correlations between ip^ and ipk+j are given by the j th power of the 
transition matrix T yielding 



,in(ipk+j-4>k) 



) = a(n) tr ■ T J ■ b( 



n 



(A3) 



with vectors 



b tr (n) = ( X (s) (0),...,e 



-(L-l)inAi/> (s) 



*W((L-1)AV0) 



Upon inserting Eq. (|A3|) into Eq. (|A^), we observe that the expression for h M depends 
only on j but not on k so that we may perform the outer summation to obtain 



M-l 



h M {nO) = a tr (n) 
Diagonalizing T leads to 



M 



b(n) . 



3=1 



h M (nn) = a(n) tr ■ S (M) ■ b(n) . 

Here, the diagonal matrix S^ M ) is given by 

M — 1 

C(M) 



(A4) 



for m = 1 



2 

A m 1 A m (A — 1) 

+ — 7TT7 rTTT" for 771 > 1 



(A5) 



with 



1 - A m M (A^-l) 2 
T = C L C -1 , L = diag(l > |A 2 | > . . . > |A L | 



a(n) = C tr a(n) , b(n) = C _1 b(n) . 
Inserting Eq. (|A4[) into Eq. (|Al|), we have 



5 T (nfi) = — 7— r [1 + 2 Re (a tr (n)S (M) b(n))l 
7r (r) 



(A6) 



Finally, we split the matrix S^ M -* into the parts pertaining to the discrete and the con- 
tinuous parts of the spectrum and define the functions A and B 

S (M) = diag(^i, 0, . . . , 0) + diag(0, S™,..., S { L f) , 



A(n, M) = 2 Re a tr (n) diag(0, S™, S%>) b(n) 
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B(n) = Re[ai(n)6i(n)] 



Rewriting Eq. (|A6|) accordingly, we arrive at the desired expression for the power spectral 
density 



s To (nn) 



7r(r) L 



1 + A(n, M) + (M — l)B(n) 



To see that A is bounded as M — * oo, note that A depends on M only through the 
diagonal entries of with m > 1 and |A m | < 1. For these we have 



lim | 1 



M^oo 



1 - Xr 



< oo , m > 1 . 
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FIG. 1. (a) Conditional ISI distributions for <fi = — 7r/6 (solid) and <p = ir/6 (dashed); other 
parameters ji = 0.95, q = 0.048, = 0.05-7T, D = 6 ■ 10~ 5 . T = 40 is the stimulus period. The 
refractory mode at small r is present only for 4> = —tt/6. The small modes around 2T correspond 
to the probability of skipping a period, (b) The same distributions as in (a), but now plotted vs. 
phase, if> = [Qt + 4>] mod 2ir, shifted to [— tt, it]. The modes at w — 7r/25 coincide for the first and 
second stimulus period, while the refractory mode is clearly set apart. 
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FIG. 2. (a) Graphic representation of the Markov chain iteration given by Eq. (U). The 
dashed line is the matrix diagonal. Probability is given by grayscale as indicated by the colorbar. 
(b) Evolution of an initially uniform phase distribution under subsequent multiplications with T, 
from right to left. See text for details. Stimulus parameters: \i = 0.95, q = 0.05, ft = 0.02-7T, 
D = 1.3- 10~ 4 . 
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FIG. 3. Phase transition matrices T (a, c, e) and corresponding stationary phase distributions 
(b, d, f) for stimulus frequency f2 = 0.05-7T at three different noise intensities D = 6.2 • 10~ 6 
(a, b), D = 7.0 ■ 1(T 5 (c, d), and D = 4.8 • 1(T 3 (e, f); other parameters: \i = 0.95, q = 0.05. 
The grayscale is the same for all matrices, white indicating vanishing probability. Error bars in 
the phase distributions indicate standard error of mean from simulated trains of 20,000 spikes. 
Observe the different scalings of the ordinate. 
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0.06 




FIG. 4. Stationary phase distributions f° r stimulus frequency f2 = 0.57T and noise intensi- 
ties D = 7.8- 10" 4 (solid), D = 4.8 • 10~ 3 (dashed), and D = 3.0 • 1(T 2 (dash-dotted). Everything 
else is as in Fig. ||. 



0.3 - 
0.25 - 



0.2 



0.15 



0.1 



0.05 



^» J. 




20 



40 



60 



80 



FIG. 5. Stationary ISI distributions for slow stimuli. Noise intensities are D = 6.2 -10 -6 (solid), 
D = 7.0 • 10~ 5 (dashed), and D = 4.8 • 10 -3 (dash-dotted). All other parameters are as in Fig. ||, 
error bars again indicate simulation results. 



20 



0.4 



~i 1 1 1 1 1 r 




FIG. 6. Stationary ISI distributions for fast stimuli. Noise intensities are D = 7.8- 10~ 4 (solid), 
D = 4.8 • 10~ 3 (dashed), and D = 3.0 • 10 -2 (dash-dotted). All parameters are as in Fig. ||, and 
error bars are from simulations. 
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FIG. 7. Power spectral density from an observation time of T a = 200 for the same stimulus 
as in Fig. ||(b, e). Circles indicate results at stimulus harmonics from the Markov chain analysis, 
while the drawn out line is obtained by FFT from a simulated train of 20,000 spikes. Ticks on the 
abscissa mark multiples of the stimulus period O = 0.05-7T. 
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FIG. 8. Signal-to-noise ratio vs. noise amplitude for three different stimulus frequencies: 
Q = 0.1-7T (dashed), Q = 0.33-7T (solid) and 17 = 0.5-7T (dash-dotted). Other parameters are fx = 0.95 
and q = 0.05. Error bars show s.e.m. from simulated trains of 20,000 spikes. 
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FIG. 9. Inter-spike-interval distributions for the resonance frequency f2 r and three noise inten- 
sities D = 1.3 • 10~ 4 (solid), D = 7.8 • 10~ 4 (dashed) and D = 4.8 • 10~ 3 (dash-dotted); other 
parameters [i = 0.95, q = 0.05. 
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FIG. 10. Signal-to-noise ratio from the phenomenological model of Eq. (^). All parameters 
are as in Fig. || with stimulus frequencies f2 = 0.17T (dashed), O, = 0.337T (solid) and f2 = 0.57T 
(dash-dotted). 
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